#!/usr/bin/perl -w
use strict;
use Bio::AlignIO; # multiple alignment parser and writer
use Bio::DB::Fasta; # database indexing and retrieval of sequences
use Getopt::Long; # cmdline arguments
use Bio::Align::Utilities qw(:all);
use Bio::Seq;

# input 
# 1. multi-fasta multiple alignment (FASTA with gaps) in protein-space (protein alignment)
# 2. Coding Sequence database (CDS)

my ($cds_db,$input,$output);

my $aln_format = 'fasta';
GetOptions(
	   'db:s'       => \$cds_db,
	   'o|output:s' => \$output,
	   'i|input:s'  => \$input, # multiple alignment input
	   
	   );


my $db = Bio::DB::Fasta->new($cds_db);
my $in = Bio::AlignIO->new(-format => $aln_format,
			   -file   => $input);
my $aa_aln = $in->next_aln;
if( ! $aa_aln ) {
    # sanity check
    die("couldn't parse an alignment out of $input\n");
}

my %dnaseqs;
for my $seq ( $aa_aln->each_seq ) { # walk through each sequence in the alignment (each domain instance)
    my $id = $seq->id; # sequence ID of the form BDET_1234.dom1_3
    my ($geneid,$domaininfo) = split(/\./,$id,2);
    my ($protein_coords) = split(/\s+/,$seq->description);
    my ($aa_start, $aa_stop) = split(/\.\./,$protein_coords);
    my $cds_start = ($aa_start - 1) * 3 + 1;
    my $cds_stop  = ($aa_stop - 1) * 3 + 1;
    my $cdsseq = $db->seq($geneid, $cds_start, $cds_stop);
    my $dnaseqobj = Bio::Seq->new(-id   => $id,
				  -desc => sprintf("%d..%d ",$cds_start,$cds_stop, $seq->description),
				  -seq  => $cdsseq,
				  );    
    $dnaseqs{$id} = $dnaseqobj;
}

# this step projects the alignment of the proteins onto the CDS sequence fragments
my $dna_aln = &aa_to_dna_aln($aa_aln,\%dnaseqs);
my $out;
if( $output ) {
    $out = Bio::AlignIO->new(-format => $aln_format,
			     -file   => ">$output");
} else {
    $out = Bio::AlignIO->new(-format => $aln_format);
}

$out->write_aln($dna_aln);
